MODULE SRFSN_LWIMPMLS_MOD
CONTAINS
SUBROUTINE SRFSN_LWIMPMLS(KIDIA  ,KFDIA  ,KLON   ,PTMST,        &
 & PSSNM1M ,PTSNM1M ,PRSNM1M ,PSSNTOTM1M, PRSNTOTM1M,PTSAM1M,PHLICEM1M,              &
 & PSLRFL  ,PSSRFLTI,PFRTI   ,PAHFSTI ,PEVAPTI,               &
 & PSSFC   ,PSSFL   ,PEVAPSNW,                                &
 & PTSFC   ,PTSFL   ,                                         &
 & PAPRS   ,PWSAM1M ,                                         &
 & PSSDP3,&
 & YDCST   ,YDVEG   ,YDSOIL  ,YDFLAKE ,                       &
 & PTSN    ,PGSN )


USE PARKIND1  , ONLY : JPIM, JPRB
USE YOMHOOK   , ONLY : LHOOK, DR_HOOK, JPHOOK
USE YOS_CST   , ONLY : TCST
USE YOS_VEG   , ONLY : TVEG
USE YOS_SOIL  , ONLY : TSOIL
USE YOS_FLAKE , ONLY : TFLAKE
USE YOMSURF_SSDP_MOD

#ifdef DOC
! (C) Copyright 2011- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.
!**** *SRFSNS_LWIMP* - CONTAINS SNOW PARAMETRIZATION
!
!     PURPOSE.
!     --------
!          COMPUTES CHANGES IN SNOW TEMPERATURE

!**   INTERFACE.
!     ----------
!          *SRFSN_SLWIMP* IS CALLED FROM *SURFTSTPS*.

!     PARAMETER   DESCRIPTION                                    UNITS
!     ---------   -----------                                    -----

!     INPUT PARAMETERS (INTEGER):
!    *KIDIA*      START POINT
!    *KFDIA*      END POINT
!    *KLON*       NUMBER OF GRID POINTS PER PACKET

!     INPUT PARAMETERS (REAL):
!    *PTMST*      TIME STEP                                      S

!     INPUT PARAMETERS AT T-1 OR CONSTANT IN TIME (REAL):
!    *PSSNM1M*    SNOW MASS (per unit area)                    kg/m**2/s
!    *PTSNM1M*    SNOW TEMPERATURE                               K
!    *PRSNM1M*    SNOW DENSITY                                 KG/M3
!    *PTSAM1M*    SOIL TEMPERATURE                               K
!    *PHLICEM1M*  LAKE ICE THICKNESS                             m
!    *PSSRFLTI*   NET SHORTWAVE RADIATION AT THE SURFACE,
!                  FOR EACH TILE                                 W/M**2
!    *PSLRFL*     NET LONGWAVE  RADIATION AT THE SURFACE         W/M**2
!    *PFRTI*      TILE FRACTIONS                                 -
!    *PAHFSTI*    TILE SENSIBLE HEAT FLUX                      W/M**2
!    *PEVAPTI*    TILE EVAPORATION                             KG/M**2/S
!    *PSSFC*      CONVECTIVE  SNOW FLUX AT THE SURFACE         KG/M**2/S
!    *PSSFL*      LARGE SCALE SNOW FLUX AT THE SURFACE         KG/M**2/S
!    *PEVAPSNW*   EVAPORATION FROM SNOW UNDER FOREST           KG/M2/S
!    *PTSFC*      Convective Throughfall at the surface        KG/M**2/S
!    *PTSFL*      Large Scale Throughfall at the surface       KG/M**2/S

!     PARAMETERS AT T+1 :
!    *PTSN*       SNOW TEMPERATURE                               K

!    FLUXES FROM SNOW SCHEME:
!    *PGSN*       GROUND HEAT FLUX FROM SNOW DECK TO SOIL     W/M**2   (#)

! (#) THOSE TWO QUANTITIES REPRESENT THE WHOLE GRID-BOX. IN RELATION
!       TO THE DOCUMENTATION, THEY ARE PGSN=Fr_s*G_s

!     METHOD.
!     -------
!     Based on the original snow (as in ERA-40) with the following updates:
!     - Liquid water as a diagnostics
!     - Interception of rainfall

!     EXTERNALS.
!     ----------

!     REFERENCE.
!     ----------
!          SEE SOIL PROCESSES' PART OF THE MODEL'S DOCUMENTATION FOR
!     DETAILS ABOUT THE MATHEMATICS OF THIS ROUTINE.

!     Original:
!     ---------
!          Simplified version based on SRFSN_LWIMP
!     M. Janiskova              E.C.M.W.F.     26-07-2011  

!     Modifications
!     -------------
!     I. Ayan-Miguez (BSC) Sep 2023 Added PSSDP3 object for spatially distributed parameters 
!     ------------------------------------------------------------------
#endif

IMPLICIT NONE

! Declaration of arguments

INTEGER(KIND=JPIM), INTENT(IN)   :: KIDIA
INTEGER(KIND=JPIM), INTENT(IN)   :: KFDIA
INTEGER(KIND=JPIM), INTENT(IN)   :: KLON

REAL(KIND=JPRB),    INTENT(IN)   :: PTMST
REAL(KIND=JPRB),    INTENT(IN)   :: PSSNM1M(:)
REAL(KIND=JPRB),    INTENT(IN)   :: PTSNM1M(:)
REAL(KIND=JPRB),    INTENT(IN)   :: PRSNM1M(:)
REAL(KIND=JPRB),    INTENT(IN)   :: PSSNTOTM1M(:)
REAL(KIND=JPRB),    INTENT(IN)   :: PRSNTOTM1M(:)
REAL(KIND=JPRB),    INTENT(IN)   :: PTSAM1M(:,:)
REAL(KIND=JPRB),    INTENT(IN)   :: PHLICEM1M(:)
REAL(KIND=JPRB),    INTENT(IN)   :: PSLRFL(:)
REAL(KIND=JPRB),    INTENT(IN)   :: PSSRFLTI(:,:)
REAL(KIND=JPRB),    INTENT(IN)   :: PFRTI(:,:)
REAL(KIND=JPRB),    INTENT(IN)   :: PAHFSTI(:,:)
REAL(KIND=JPRB),    INTENT(IN)   :: PEVAPTI(:,:)
REAL(KIND=JPRB),    INTENT(IN)   :: PSSFC(:)
REAL(KIND=JPRB),    INTENT(IN)   :: PSSFL(:)
REAL(KIND=JPRB),    INTENT(IN)   :: PEVAPSNW(:)
REAL(KIND=JPRB),    INTENT(IN)   :: PTSFC(:)
REAL(KIND=JPRB),    INTENT(IN)   :: PTSFL(:)
REAL(KIND=JPRB),    INTENT(IN)   :: PSSDP3(:,:,:)
TYPE(TCST),         INTENT(IN)   :: YDCST
TYPE(TVEG),         INTENT(IN)   :: YDVEG
TYPE(TSOIL),        INTENT(IN)   :: YDSOIL
TYPE(TFLAKE),       INTENT(IN)   :: YDFLAKE

REAL(KIND=JPRB),    INTENT(OUT)  :: PTSN(:)
REAL(KIND=JPRB),    INTENT(OUT)  :: PGSN(:)

REAL(KIND=JPRB),    INTENT(IN)   :: PAPRS(:), PWSAM1M(:,:)

!      LOCAL VARIABLES

LOGICAL            :: LLNOSNOW(KLON)

REAL(KIND=JPRB),DIMENSION(KLON) :: ZFRSN,ZIFRSN,ZSSFC,ZSSFL,ZLWC,ZDSNR,ZDSNRP,&
 &        ZTSFC,ZTSFL,ZLICE
REAL(KIND=JPRB) :: ZFUNC,ZDFUNC,ZDSN,ZDTH,ZDTHI,ZDTM,&
 &        ZEXPF,ZGSN,ZHFLUXP,ZHFLUX,ZLINA,ZMSN,ZRS,ZSNRES,ZSOILRES,ZSSTAR,&
 &        ZSSTAR1,ZT0,ZTSTAR,ZTSTARI,ZHOICE,ZTMST,ZCONSA,ZCONSAMAX,ZCONSB,&
 &        ZCONS2,ZCONSLWC,ZWFLUXRF,ZHFLUXRF1,ZITEMPAMP,&
 &        ZGRIDFRAC
REAL(KIND=JPRB) :: ZCOND1,ZCOND2,ZCOND3,ZCOND4,ZCOND5,ZCOND6
REAL(KIND=JPRB) :: ZRDSNMAX, ZRDSNRESMAX

REAL(KIND=JPRB)                   :: ZSNVCOND, ZSNCOND, ZSOILDEPTH1 
REAL(KIND=JPRB),DIMENSION(KLON)   :: ZSOILCOND

REAL(KIND=JPRB) :: ZFF, ZWU, ZLIC, ZLWT, ZLAMBDASAT, ZKERSTEN, ZINVWSAT, ZCOND

INTEGER(KIND=JPIM) :: JL

REAL(KIND=JPHOOK) :: ZHOOK_HANDLE

! Inlude functions!
#include "fcsurf.h"

!     ------------------------------------------------------------------

IF (LHOOK) CALL DR_HOOK('SRFSN_LWIMPMLS_MOD:SRFSN_LWIMPMLS',0,ZHOOK_HANDLE)
ASSOCIATE(RDAY=>YDCST%RDAY, RLMLT=>YDCST%RLMLT, RLSTT=>YDCST%RLSTT, &
 & RLVTT=>YDCST%RLVTT, RPI=>YDCST%RPI, RTT=>YDCST%RTT, &
 & LEFLAKE=>YDFLAKE%LEFLAKE, RH_ICE_MIN_FLK=>YDFLAKE%RH_ICE_MIN_FLK, &
 & RALAMSN=>YDSOIL%RALAMSN, RDSNMAX=>YDSOIL%RDSNMAX, RFRSMALL=>YDSOIL%RFRSMALL, &
 & RFRTINY=>YDSOIL%RFRTINY, RHOCI=>YDSOIL%RHOCI, RHOICE=>YDSOIL%RHOICE, &
 & RLAMICE=>YDSOIL%RLAMICE, RLWCSWEA=>YDSOIL%RLWCSWEA, &
 & RLWCSWEB=>YDSOIL%RLWCSWEB, RLWCSWEC=>YDSOIL%RLWCSWEC, RTAUA=>YDSOIL%RTAUA, &
 & RTAUF=>YDSOIL%RTAUF, RTEMPAMP=>YDSOIL%RTEMPAMP, &
 & RVLAMSK_DESERT=>YDVEG%RVLAMSK_DESERT, &
 & SNHCONDAV=>YDSOIL%SNHCONDAV, &
 & SNHCONDBV=>YDSOIL%SNHCONDBV, SNHCONDCV=>YDSOIL%SNHCONDCV, &
 & SNHCONDPOV=>YDSOIL%SNHCONDPOV, &
 & RKERST1=>YDSOIL%RKERST1, RKERST2=>YDSOIL%RKERST2, RKERST3=>YDSOIL%RKERST3, &
 & RLAMBDADRY=>YDSOIL%RLAMBDADRY, RLAMBDADRYM3D=>PSSDP3(:,:,SSDP3D_ID%NRLAMBDADRYM3D), &
 & RLAMBDAICE=>YDSOIL%RLAMBDAICE, RLAMBDAWAT=>YDSOIL%RLAMBDAWAT, &
 & RLAMSAT1=>YDSOIL%RLAMSAT1, RLAMSAT1M3D=>PSSDP3(:,:,SSDP3D_ID%NRLAMSAT1M3D), RSIMP=>YDSOIL%RSIMP, &
 & RWSAT=>YDSOIL%RWSAT, RWSATM3D=>PSSDP3(:,:,SSDP3D_ID%NRWSATM3D))

!*         1.    SET UP SOME CONSTANTS.
!                --- -- ---- ----------
!*               PHYSICAL CONSTANTS.
!                -------- ----------
! RESISTANCE FOR HEAT CONDUCTION IN HALF TOP SOIL LAYER
!  (ESTIMATE FROM SKIN LAYER CONDUCTIVITY FOR BARE SOIL; W/M2K)


ZHOICE=1.0_JPRB/RHOICE
ZSOILRES=1.0_JPRB/RVLAMSK_DESERT
ZT0=RTT
ZTMST=1.0_JPRB/PTMST
ZEXPF=EXP(-RTAUF*PTMST/RDAY)
ZLINA=RTAUA*PTMST/RDAY
ZITEMPAMP=1.0_JPRB/RTEMPAMP

ZRDSNMAX=0.05_JPRB
ZRDSNRESMAX=100._JPRB
ZSOILDEPTH1=YDSOIL%RDAW(1)

!*   FIND LAKE POINTS WITH ICE COVER
DO JL=KIDIA,KFDIA
  IF ( PHLICEM1M(JL) > RH_ICE_MIN_FLK  .AND. LEFLAKE ) THEN
    ZLICE(JL)=1._JPRB
  ELSE
    ZLICE(JL)=0._JPRB
  ENDIF
ENDDO

!     ------------------------------------------------------------------
!*         2. NEW SNOW T AND MASS INCLUDING GROUND HEAT FLUX AND MELTING.
!             -----------------------------------------------------------

DO JL=KIDIA,KFDIA
  ZFRSN(JL)=MAX(PFRTI(JL,5)+PFRTI(JL,7),RFRTINY)
  ZIFRSN(JL)=1.0_JPRB/ZFRSN(JL)
  IF (ZFRSN(JL) < RFRSMALL) THEN
    LLNOSNOW(JL)=.TRUE.
  ELSE
    LLNOSNOW(JL)=.FALSE.
  ENDIF
  ZGRIDFRAC=(PFRTI(JL,3)+PFRTI(JL,4)+PFRTI(JL,5)&
   & +PFRTI(JL,6)+PFRTI(JL,7)+PFRTI(JL,8))

  IF ( LEFLAKE ) THEN
    IF ( PFRTI(JL,9) .EQ. 1._JPRB ) THEN
      ZGRIDFRAC=0._JPRB ! FULL COVER OF LAKE --NO SNOW THERE !
    ELSE
      ZGRIDFRAC=ZGRIDFRAC+PFRTI(JL,9)
    ENDIF
  ENDIF

  ZSSFC(JL)=ZGRIDFRAC*PSSFC(JL)
  ZSSFL(JL)=ZGRIDFRAC*PSSFL(JL)
  ZTSFC(JL)=ZGRIDFRAC*PTSFC(JL)
  ZTSFL(JL)=ZGRIDFRAC*PTSFL(JL)

ENDDO

DO JL=KIDIA,KFDIA
  IF (LLNOSNOW(JL)) THEN
    ZSSTAR=PSSNM1M(JL)+PTMST*(ZSSFC(JL)+ZSSFL(JL))
    PTSN(JL)=ZT0
    PGSN(JL)=0.0_JPRB
    ZLWC(JL)=0.0_JPRB
  ELSE
    ZFF=0.0_JPRB
    ZWU=PWSAM1M(JL,1)*(1.0_JPRB-ZFF)
    ZLWT=RLAMBDAWAT**ZWU
    ZINVWSAT=1.0_JPRB/RWSATM3D(JL,1)
    ZLIC=RLAMBDAICE**(RWSATM3D(JL,1)-ZWU)
    ZLAMBDASAT=RLAMSAT1M3D(JL,1)*ZLIC*ZLWT
    ZCOND=MAX(RKERST1,PWSAM1M(JL,1)*ZINVWSAT)
    ZKERSTEN=RKERST2*LOG10(ZCOND)+RKERST3
    ZSOILCOND(JL)=RLAMBDADRYM3D(JL,1)+ZKERSTEN*(ZLAMBDASAT-RLAMBDADRYM3D(JL,1))


!           NET HEAT FLUX AT SNOW SURFACE; THIS IS MULTIPLIED
!           BY FRACTION BECAUSE EQUATIONS APPLY TO TOTAL SNOW MASS
!           IN THE GRID SQUARE

    ZHFLUXP=PFRTI(JL,5)*(PAHFSTI(JL,5)+RLSTT*PEVAPTI(JL,5))&
     & +PFRTI(JL,7)*(PAHFSTI(JL,7)+RLVTT*(PEVAPTI(JL,7)-PEVAPSNW(JL))&
     & +RLSTT*PEVAPSNW(JL))&
     & +PFRTI(JL,5)*PSSRFLTI(JL,5)&
     & +PFRTI(JL,7)*PSSRFLTI(JL,7)&
     & +(PFRTI(JL,5)+PFRTI(JL,7))*PSLRFL(JL)
    ZHFLUX=ZIFRSN(JL)*ZHFLUXP

! PRELIMINARY SNOW MASS
    ZSSTAR=PSSNM1M(JL)+PTMST&
     & *(ZSSFC(JL)+ZSSFL(JL)+PFRTI(JL,5)*PEVAPTI(JL,5)&
     & +PFRTI(JL,7)*PEVAPSNW(JL))
    ! CORRECT FOR NEGATIVE VALUES OF SNOW MASS AFTER P-E
    IF (ZSSTAR < 0.0_JPRB) THEN
!      ZSSTAR=MAX(ZSSTAR,0.0_JPRB)
      ZSSTAR=0.0_JPRB
    ENDIF

! SNOW LIQUID WATER - IMPLICIT
    ! SNOW DEPTH (REAL)
    ZDSNR(JL)=PSSNM1M(JL)/(PRSNM1M(JL)*ZFRSN(JL))
    ! SNOW LIQUID WATER CAPACITY
    ZCONSLWC=MAX(0._JPRB,(RLWCSWEA-PRSNM1M(JL))/RLWCSWEA)
    ZLWC(JL)=PSSNM1M(JL)*(RLWCSWEB+(RLWCSWEC-RLWCSWEB)*ZCONSLWC)
    ! LWC IS RESCALLED CONSIDERING EFFECTIVE DEPTH (ZRDSNMAX)
    ! ZLWC(JL)=ZLWC(JL)*MIN(ZDSNR(JL),ZRDSNMAX)/ZDSNR(JL)
    IF (ZDSNR(JL) < ZRDSNMAX) THEN
      ZDSNRP(JL) = ZDSNR(JL)
    ELSE
      ZDSNRP(JL) = ZRDSNMAX
    ENDIF
    ZLWC(JL)=ZLWC(JL)*ZDSNRP(JL)/ZDSNR(JL)
    ! ANALYTICAL FUNCTIONS - SNOW LIQUID WATER
    IF ( PTSNM1M(JL) < ZT0-0.5_JPRB*RTEMPAMP ) THEN
      ZFUNC=0._JPRB
      ZDFUNC=0._JPRB
    ELSEIF(PTSNM1M(JL) < ZT0) THEN
      ZFUNC=1._JPRB+SIN(RPI*(PTSNM1M(JL)-ZT0)*ZITEMPAMP)
      ZDFUNC=COS(RPI*(PTSNM1M(JL)-ZT0)*ZITEMPAMP)*RPI*ZITEMPAMP
    ELSE
      ZFUNC=0._JPRB
      ZDFUNC=RPI*ZITEMPAMP
    ENDIF

! TERMS FOR IMPLICIT TEMPERATURE EQUATION
    ZCONSB=RLMLT*ZLWC(JL)*ZDFUNC
    ZCONSAMAX=(PRSNM1M(JL)*RHOCI*ZRDSNMAX)*ZHOICE
    ZCOND1=(PSSNM1M(JL)*RHOCI*ZHOICE)*ZIFRSN(JL)
    ZCONSA=MIN(ZCOND1,ZCONSAMAX)
    ZCONS2=PTMST/(ZCONSA+ZCONSB)

! REAL SNOW DEPTH
    ZCOND2=PSSNTOTM1M(JL)/(PRSNTOTM1M(JL)*ZFRSN(JL))
    ZDSN=MIN(ZCOND2,ZRDSNRESMAX)
    ! RESISTANCE FOR HEAT FLUX BETWEEN SNOW AND SOIL LAYER
!   ZSNRES=0.5_JPRB*ZDSN/(RLAMICE*(PRSNM1M(JL)*ZHOICE)**RALAMSN)
!   ZRS=1.0_JPRB/(ZSNRES+ZSOILRES)
    ZSNVCOND=(SNHCONDPOV/PAPRS(JL))*MAX(0._JPRB,(SNHCONDAV-SNHCONDBV/(PTSNM1M(JL)-SNHCONDCV)))  
    ZSNCOND = FSNTCOND(PRSNTOTM1M(JL))+ZSNVCOND 
    ZSNRES  = 0.5_JPRB*ZDSN/ZSNCOND
    ZSOILRES= 1.0_JPRB*ZSOILDEPTH1/ZSOILCOND(JL)

    ZRS=1.0_JPRB/(ZSNRES+ZSOILRES)

!  ACCOUNT FOR THROUGHFALL INTERCEPTION
    ! mass advected
    ZWFLUXRF=(PFRTI(JL,5)+PFRTI(JL,7))*(ZTSFC(JL)+ZTSFL(JL))
    ! phase changes
    ZCOND3=ZWFLUXRF*PTMST
    ZCOND4=ZLWC(JL)*(1.0_JPRB-ZFUNC)
    ZWFLUXRF=MIN(ZCOND3,ZCOND4)
    ZCOND5=ZCONSA*(ZT0-PTSNM1M(JL))*ZTMST+ZT0*ZRS
    ZCOND6=RLMLT*ZWFLUXRF*ZTMST
    ZHFLUXRF1=MIN(ZCOND5,ZCOND6)
    ZHFLUXRF1=MAX(ZHFLUXRF1,0.0_JPRB)
    !mass intercepted - recalculate from ZHFLUXRF1
    ZWFLUXRF=ZHFLUXRF1/RLMLT
    ZSSTAR=ZSSTAR+ZWFLUXRF*PTMST
! SOLVE IMPLICIT TEMPERATURE EQUATION
    ZTSTARI=1.0_JPRB/(1.0_JPRB+ZCONS2*ZRS)
    ZTSTAR=(PTSNM1M(JL)+ZCONS2*(ZHFLUX+ZHFLUXRF1+PTSAM1M(JL,1)*ZRS)) &
     & *ZTSTARI
    IF (ZTSTAR <= ZT0) THEN
! NO MELTING OF SNOW
      PTSN(JL)=ZTSTAR
      PGSN(JL)=(ZTSTAR-PTSAM1M(JL,1))*ZRS
    ELSE
! MELTING OF SNOW:
! HEATING PART OF THE TIME STEP
      ZDTHI=1.0_JPRB/(ZHFLUX+ZHFLUXRF1-(ZT0-PTSAM1M(JL,1))*ZRS)
      ZDTH=((ZT0-PTSNM1M(JL))*(ZCONSA+ZCONSB))*ZDTHI 
! MELTING PART OF THE TIME STEP
      ZDTM=PTMST-ZDTH
      ZGSN=(ZT0-PTSAM1M(JL,1))*ZRS
      ZMSN=(ZHFLUX+ZHFLUXRF1-ZGSN)/RLMLT
      ZSSTAR1=ZSSTAR-ZDTM*ZMSN*ZFRSN(JL)
      IF (ZSSTAR1 >= 0.0_JPRB) THEN
! NOT ALL THE SNOW HAS BEEN MELTED
        PTSN(JL)=ZT0
        PGSN(JL)=(ZT0-PTSAM1M(JL,1))*ZRS
      ELSE
! ALL THE SNOW HAS BEEN MELTED; COMPUTE TIME IT TOOK TO MELT
        PTSN(JL)=ZT0
        PGSN(JL)=( ((ZT0-PTSAM1M(JL,1))*ZRS)*(ZDTM+ZDTH)&
         & +(ZHFLUX+ZHFLUXRF1)*(PTMST-ZDTM-ZDTH) )*ZTMST
      ENDIF
    ENDIF

  ENDIF
ENDDO

!*         5. NORMALIZE QUANTITIES TO THE GRID-SQUARE.
!             ----------------------------------------

DO JL=KIDIA,KFDIA
  PGSN(JL)=ZFRSN(JL)*PGSN(JL)
ENDDO

END ASSOCIATE
IF (LHOOK) CALL DR_HOOK('SRFSN_LWINPMLS_MOD:SRFSN_LWIMPMLS',1,ZHOOK_HANDLE)

END SUBROUTINE SRFSN_LWIMPMLS
END MODULE SRFSN_LWIMPMLS_MOD
